Porting a Kalman filter based track fit to NVIDIA CUDA

Finding particle trajectories is usually the most time consuming part of modern experiments in high energy physics. In many present experiments with high track densities and complicated event topologies a Kalman filter based track fit is used already at this combinatorial part of the event reconstruction. Therefore speed of the track fitting algorithm becomes very important for the total processing time. In 2007 a Kalman filter based track fitting algorithm of the CBM experiment [1] had been ported to SSE and the Cell SPE [2]. On the CPU a speedup of 10000 compared to the initial version was achieved. To gain this speedup two major changes to the existing Kalman filter implementation were done. The magnetic field map was replaced by a functional approximation based on fourth order polynomials. In addition the algorithm was tuned to work with single precision calculations only, without the need for compensating computations or numerical instabilities. Today CPUs are no longer able to increase their peak performance by increasing clock speeds. Modern GPUs however continue to increase their peak performance utilising manycore architectures, with NVIDIA GT200 based GPUs nearly reaching 1 TFlop/s in single precision [3]. Thus they are a promising candidate for further acceleration of the algorithm. This performance comes at a cost. Different than on CPUs access to the main memory is not cached, all threads running on one multiprocessor, which is similar to a core of a GPU, share a small 16 kB processor local storage called shared memory which needs to be explicitly programmed. To compensate for that each multiprocessor has a large register file of 16384 registers and can have 1024 threads concurrently active. On eight ALUs it processes warps of 32 threads in 4 cycles, switching between warps without costs, allowing to hide memory accesses by calculations. A GT200 chip contains 30 of these multiprocessors. In contrast to the single instruction multiple data (SIMD) based architectures of the NVIDIA GPUs are based on a single instruction multiple threads (SIMT) model, which means each ALU is connected with an own instruction counter. Therefore the algorithm needs to be parallized on a thread level instead of the below as in the SIMD case. While in the SIMD case we have one thread that works on a vector of numbers, in SIMT we have multiple threads working on scalars. However in both cases a number of arithmetic units is fed by the same instruction decoder, making it important for the code not to rely on different code branches depending on the data of each thread. The original port relied on operator overloading to have only one implementation of the algorithm for multiple plat-